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Abstract 

Random Boolean networks (RBNs) are frequently employed for modelling complex systems 
driven by information processing, e.g. for gene regulatory networks (GRNs). Here we propose 
a hierarchical adaptive RBN (HARBN) as a system consisting of distinct adaptive RBNs - 
subnetworks - connected by a set of permanent interlinks. Information measures and internal 
subnetworks topology of HARBN coevolve and reach steady-states that are specific for a given 
network structure. We investigate mean node information, mean edge information as well as 
a mean node degree as functions of model parameters and demonstrate HARBNs ability to 
describe complex hierarchical systems. 

Keywords: coevolution, information processing, topology, adaptive random Boolean networks, hierar¬ 
chical 


1 Introduction 

Complex networks are frequently used to describe evolving, constantly changing objects that 
consist of elements and complicated functional relations between them. Simplified models were 
created to study and compare certain network properties. One of the examples are random 
Boolean networks (RBNs). RBNs are generic [T! and that is why they have been applied in 
many different fields. 

Gene regulatory networks (GRNs) are models describing the structure and behavior of the 
transcriptional network responsible for regulating the gene expression in a cell [2, 3]. GRNs 
are both robust and stochastic. GRNs must respond to diverse stimuli and they are highly 
modular, e.g. hierarchical regulatory interactions in transcriptional networks in yeast. Analysis 
of GRNs provided evidence that diverse stimuli may cause modifications of gene interactions 
and network topology. During GRN evolution edges are created and destroyed, and in this way 
the gene expression is altered. 
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In [1] Liu and Bassler proposed the RBN modification, where network topology is time 
dependent and the observed coevolution encompasses some effects noticed in GRN. Here, we 
extend the original Liu and Bassler’s model, most of all by incorporating the concept of hier¬ 
archy. In our approach there are no robust genes, but a part of network edges are assumed to 
be robust, i.e. they stay unchanged during the system evolution. 

The paper is organised as follows: in section [2] we describe our model, i.a. we present 
the algorithm of network evolution (section 12.21) . Section [3] presents the applied information 
measures. Section [|] contains simulations description, results and discussion. Finally, in section 
[5] we conclude our studies. 

2 Model of Hierarchical Adaptive Random Boolean Net¬ 
work (HARBN) 

Let us consider a system of coupled nodes in a form of a directed network that corresponds to a 
GRN. For simplicity we will assume that internal nodes’ variables a n (t) are 0 or 1 and they can 
evolve in time according to randomly chosen boolean functions /). Arguments of the function 
f n are internal variables <r m ( t ) of all nodes m such that there is a connection from a node m to 
the node n. 

During simulations of adaptive RBN (ARBN) changes of network state and topology are 
measured. The simulation consists of defined a priori number of epochs. In each epoch the 
network’s attractor is found (i.e a periodic orbit of variables cr n (t) that is reached by the system 
after a transient time), one node is randomly chosen and according to the Activity-dependent 
Rewiring Rule (ADRR) (described in section 12.11) a number of incoming connections to this 
node is changed. 

We extended Liu and Bassler’s model. First, we consider a small number of ARBNs, called 
subnetworks. Different subnetworks are sparsely connected by directed edges, called interlinks. 
We build our model in such a way that the considered number of interlinks is smaller than 
the number of edges within the subnetworks. Therefore, the created network is hierarchical 
in terms of connection density, because certain groups of nodes are much denser connected in 
internal blocks (subnetworks) as compared to the rest of the network. During a single epoch the 
network’s attractor is found and for each subnetwork a node is chosen randomly. Connections 
of such a node are changed according to the ADRR. The edges may arise only between nodes in 
the same subnetwork. The interlinks are permanent edges: they are created at the beginning of 
each simulation and they cannot be changed. Secondly, we also alter the ADRR by introducing 
a resilience parameter a (see section 12.11) . In the original ADRR the resilience parameter did 
not exist and it corresponds to 0 in our approach. 


2.1 Activity-Dependent Rewiring Rule 

Let us define the network resilience parameter that will decide on the type of changes in network 
topology as a where a G [0,0.5]. 

If the node’s mean state ( a n ) during the attractor is not higher than a or at least (1 — a), 
then this node is considered to be frozen and one new incoming edge is added to this node. The 
new edge starts in a random node chosen from the group of nodes from the same subnetwork 
that does not already have a link to the considered node. Otherwise, this node is considered to 
be active and one of its edges is randomly chosen and deleted. 


2 


Coevolution of Information Processing and Topology 


Gorski, Czaplicka and Holyst 


For non-zero values of parameter a the original ADRR is extended so as nodes with a mean 
state higher than 0 or smaller than 1 are considered to be frozen. 

2.2 Adaptive Algorithm for Hierarchical RBNs 

Our algorithm for the coevolution of hierarchical RBNs is as follows: 

1. Generate M uniform Boolean networks (subnetworks) containing Nm nodes each with Ki 
directed incoming edges starting in a randomly chosen group of nodes belonging to the 
same subnetwork. Then generate Km edges between subnetworks (interlinks). (For each 
node generate necessary Boolean functions with a bias parameter p = 1/2.) 

2 . Generate a random initial state S(0) of all internal node variables a n { 0) and find the 
network’s attractor length using the following algorithm [3]: 

a. Define comparison moments. Here: 

T = {0; 100; 1,000; 10,000; 100,000}. Set k = 0 and i = 1. 

b. Synchronically update the states of the network’s nodes: S(i). If S(i) equals S(Tk) the 

attractor length is (i — T*,). Go to point 3. Otherwise continue. 

c. If i equals T/c+i, increment k. If k is higher than 4 end this subalgorithm with the 

attractor length equal to (T 4 — T 3 ). Otherwise, increment i and go back to point 2b. 

3. Choose a node from each subnetwork and calculate its mean state during one attractor 
cycle. 

4. According to the ADRR, change the topology of each subnetwork. Do not delete in¬ 
terlinks. In case there are no connections to delete, randomly choose another node and 
repeat this point. 

5. Generate new Boolean functions for each node. 

6 . If the predefined number of epochs is not reached, go back to point 2. 

Remarks: 

• The interlinks are generated as follows: first the highest possible number of interlinks is 
equally distributed between ordered pairs of subnetworks (due to directed edges each pair 
of subnetworks appear twice); secondly, the remaining interlinks are generated between 
random pairs. 

• The described method allows finding the attractor length equal to 90,000 at most. 

• Advancing through the steps 2-6 is called an epoch (after [4]). Updating the states of all 
nodes in step 2 b. is called an iteration. 

3 Information Measures of RBNs 

Measuring information amount transmitted by a network is an important tool that facilitates 
exploring system’s features. There are many approaches to do it, see e.g. mi- Here we define 
network activity information I as a sum of transformed activities of each node. Activity A n of 
a node n is a number of changes of the node’s state during the network attractor divided by 
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Figure 1: Example structure after 1000 epochs of a network 4x15x10. 


the attractor length. There is no information stored in frozen nodes (nodes with A n = 0) and 
in nodes whose state changes in every iteration ( A n = 1). Such nodes should not contribute to 
the network information. On the other hand, nodes which change their state the same number 
of times as they keep the old state ( A n = 0.5) contribute mostly because of their unexpected 
behavior. Node activity can be identified with probability of changing the node’s state. Then 
the activity A n links with the parameter a described in [5] as system’s ’’self-overlap”: A n = 1—a. 
Here we introduce a natural definition of node activity information as: 

In = A n log(A n ) - (1 - A n ) log(l - A n ) (1) 

The total network activity information is given as I = For brevity we further write 

network information instead of network activity information. 


4 Simulations 

The simulations of different structures of HARBN model consisted of 1000 epochs repeated 100 
times. We let MxNmxKm denote a network consisting of M subnetworks each of Nm nodes 
and linked by Km interlinks. Here we show the results for 60— and 80—node networks. Our 
simulations were conducted in 3 parts. In the beginning we reproduced the results achieved in 
H (M = 1, a = 0). Then, we divided networks into 2 — 4 subnetworks and linked them by 
0 — 80 interlinks. We compared the results for ARBNs and HARBNs. At the end the system 
properties for non-zero resilience values a were explored. An example structure of simulated 
network after 1000 epochs is shown in Figure [Q 

Each realization of ARBNs and HARBNs demonstrates different structure and information 
parameters. However identical networks tend to oscillate (after the initial transient period) near 
the same mean steady-state (m.s.s.) levels. On the other hand, networks of different number 
of nodes, subnetworks or interlinks tend to reach different m.s.s. values. Therefore, for each 
network type we define: m.s.s. incoming connectivity K ss , m.s.s. network information J, m.s.s. 
node information IPN (calculated as information / per node), m.s.s. edge information I PE 
(calculated as information I per edge) and geometric m.s.s. attractor length T. If not other 
mentioned we assume an arithmetic mean. In order to determine all above parameters 200 
beginning epochs of each realization were discarded as transient periods. K ss is calculated as 
follows: 

• calculate mean incoming connectivity in an epoch —> Ki n , 

• calculate mean A, n starting from 201th epoch in a realization —> (A, n ), 
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Figure 2: The mean steady-state connectivity K ss as a function of the ratio of the interlinks to 
the total number of edges for various networks. 


• calculate mean ( Ki n ) over all realizations —> K ss . 

The information parameters / and IPN are calculated likewise. I PE is computed by dividing 
IPN by K ss . 

Now we shall discuss the influence of network topology on system properties. We are aware 
of the fact that because of a small number of investigated networks additional large-scale sim¬ 
ulations will be necessary to confirm our findings. 

4.1 Mean Steady-State Connectivity 

Let X denote the ratio of the interlinks to the total number of edges. For small X values 
different network types differentiate themselves in terms of K ss by the size of a subnetwork 
(Figure [2| and K ss decreases with the growth of the subnetwork size. With the increase of 
X the subnetworks become more and more linked. For high X values different network types 
start to differentiate themselves by the size of the whole network, i.e. subnetworks are strongly 
mutually dependent. We estimate that by X ss 0.40 distinct subnetworks do not exist anymore. 
The resulting K ss values are independent from the subnetwork size. 

4.2 Information per Node 

Calculated information per node IPN also exhibits different ordering for small and large X 
values (Figure [3]). For X <C 0.1 networks (apart from those with one subnetwork) differentiate 
with the increase of the number of distinct subnetworks. When X > 0.1 then the information 
IPN differentiates the networks by the size of the whole network. IPN is larger in smaller 
networks. In the interval 0 < X < 0.1, especially for 80—node networks, ambiguous behaviour 
is observable. For the 4x20 network the information IPN starts at a rather low value. A small 
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Figure 3: The mean node information IPN as a function of the ratio of the interlinks to the 
total number of edges for various networks. 


increase of the ratio X causes a quick rise of IPN with a maximum for X ss 0.05. Then until 
X ss 0.1 a moderate decrease of IPN is observable. The behaviour of IPN for 2x40 is different. 
IPN starts from a high value for X = 0 and then quickly falls to a minimum in X ss 0.05. 
Afterwards a stable increase is visible. 

4.3 Information per Edge 

The behaviour of the third network parameter — IPE (Figure [4j includes some features of 
the previous two observables. For X sa 0 networks differentiate themselves first by the number 
of distinct subnetworks (less subnetworks, higher IPE), secondly by the size of an individual 
subnetwork (with one exception for networks 1 xN). On the other hand, ordering by the network 
size for higher X is very weak and it appears that IPE is equal for all systems investigated by 
the same X value. It does not depend on the total network size or the number of subnetworks. 

4.4 Geometric Mean Attractor Length 

Figure [5] shows the relationship between the geometric mean attractor length T and the ratio 
X. The lengths of the attractors grow with the increase of the total number of nodes. Moreover, 
structures of the same size with more subnetworks tend to possess longer attractors. Dividing 
a network into subnetworks creates separate unconnected parts (X ss 0). Such structures reach 
the highest T values. It is related to the attractor search subalgorithm: the attractor is found for 
the whole network and its separate parts multiply and elongate the attractor length. Interlinks 
connect the subnetworks, but interlinks are not flexible. Interlinks significantly decrease the 
attractor length leading to higher flexibility. The smallest attractors are achieved for X ranging 
from around 0.1 to around 0.25. Further increase of the interlinks ratio leads to longer lengths 
T, as interlinks stiffen the network. 
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Figure 4: The mean edge information I PE as a function of the ratio of the interlinks to the 
total number of edges for various networks. IPE has been calculated by dividing the mean 
node information IPN by the mean steady-state in-degree connectivity K ss . 



Figure 5: The geometric mean attractor length T as a function of the ratio of the interlinks to 
the total number of edges for various networks. 
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4.5 Resilience Parameter 

In order to explore non-zero values of the resilience parameter a different network structures 
with 10 and 40 interlinks were analysed. These numbers of interlinks corresponded to X ft 0.05 
and X ft 0.25 respectively. Effects of hierarchical structures and interactions between different 
subnetworks were observed. Figure [ 6 ] shows that the parameter a differentiates networks of 
different types. First of all we can distinguish number of interlinks: there are 3 levels for 0, 
10, and 40 interlinks. Their order is primarily the result of values demonstrated in Figure 2J 
Within the networks of the same number of interlinks higher mean edge information corresponds 
to networks with larger subnetworks. Comparing 3a:20a;10 and 4x202:10 networks we can see 
the next level of ordering: higher IPE is in a smaller network. We have also observed (it 
is not shown in the paper) that non-zero values of a lead to higher K ss values. Therefore 
the information per node IPN grows with the steady-state-connectivity K ss faster than the 
information per edge I PE. 



-©-1x60 
-1-2x30x10 
-#-2x30x40 
-*-3x20x10 
-a- 3x20x40 
4x15x10 
-A- 4x15x40 
4x20x10 


Figure 6 : The mean edge information IPE as a function of the resilience parameter a for several 
different network types. IPE has been calculated by dividing the mean node information IPN 
by the mean steady-state in-degree connectivity K ss . 


5 Conclusions 

The model of hierarchical adaptive random Boolean network (HARBN) has been introduced and 
numerically explored. The system consists of subnetworks connected by fixed interlinks, where 
the internal topology of the subnetworks can evolve depending on individual node activity. We 
have observed that the main natural feature of ARBNs, i.e. their adaptability is preserved in 
HARBNs that can evolve towards stable configurations. When the ratio X of interlinks to the 
total number of edges is of the order of X ft 0.4 then the interlinks efficiently re-connect the 




































Coevolution of Information Processing and Topology 


Gorski, Czaplicka and Holyst 


whole network (as if separate parts did not exist). However such a large system is less flexible 
and displays longer attractors. The shortest attractors are observable for X ranging from 0.1 
to 0.25 depending on the network structure. 

The mean node information ( IPN ) and the mean edge information (I PE) grow in HARBN 
with the increase of the ratio X as well as with the resilience parameter a and I PE tends to 
achieve the same value for all network types in case of many interlinks. Adding a new node to 
the network decreases IPN and, moreover, leads to fewer incoming connections per node. On 
the other hand adding a new interlink increases I PE values. We conclude that the introduced 
HARBNs may successfully be used to model GRNs with a modular structure and they well 
describe the processes of evolution and speciation. 
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